home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qelg.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  5.5 KB  |  234 lines
  1. /* integration/qelg.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. struct extrapolation_table
  21.   {
  22.     size_t n;
  23.     double rlist2[52];
  24.     size_t nres;
  25.     double res3la[3];
  26.   };
  27.  
  28. static void
  29.   initialise_table (struct extrapolation_table *table);
  30.  
  31. static void
  32.   append_table (struct extrapolation_table *table, double y);
  33.  
  34. static void
  35. initialise_table (struct extrapolation_table *table)
  36. {
  37.   table->n = 0;
  38.   table->nres = 0;
  39. }
  40. #ifdef JUNK
  41. static void
  42. initialise_table (struct extrapolation_table *table, double y)
  43. {
  44.   table->n = 0;
  45.   table->rlist2[0] = y;
  46.   table->nres = 0;
  47. }
  48. #endif
  49. static void
  50. append_table (struct extrapolation_table *table, double y)
  51. {
  52.   size_t n;
  53.   n = table->n;
  54.   table->rlist2[n] = y;
  55.   table->n++;
  56. }
  57.  
  58. /* static inline void
  59.    qelg (size_t * n, double epstab[], 
  60.    double * result, double * abserr, 
  61.    double res3la[], size_t * nres); */
  62.  
  63. static inline void
  64.   qelg (struct extrapolation_table *table, double *result, double *abserr);
  65.  
  66. static inline void
  67. qelg (struct extrapolation_table *table, double *result, double *abserr)
  68. {
  69.   double *epstab = table->rlist2;
  70.   double *res3la = table->res3la;
  71.   const size_t n = table->n - 1;
  72.  
  73.   const double current = epstab[n];
  74.  
  75.   double absolute = GSL_DBL_MAX;
  76.   double relative = 5 * GSL_DBL_EPSILON * fabs (current);
  77.  
  78.   const size_t newelm = n / 2;
  79.   const size_t n_orig = n;
  80.   size_t n_final = n;
  81.   size_t i;
  82.  
  83.   const size_t nres_orig = table->nres;
  84.  
  85.   *result = current;
  86.   *abserr = GSL_DBL_MAX;
  87.  
  88.   if (n < 2)
  89.     {
  90.       *result = current;
  91.       *abserr = GSL_MAX_DBL (absolute, relative);
  92.       return;
  93.     }
  94.  
  95.   epstab[n + 2] = epstab[n];
  96.   epstab[n] = GSL_DBL_MAX;
  97.  
  98.   for (i = 0; i < newelm; i++)
  99.     {
  100.       double res = epstab[n - 2 * i + 2];
  101.       double e0 = epstab[n - 2 * i - 2];
  102.       double e1 = epstab[n - 2 * i - 1];
  103.       double e2 = res;
  104.  
  105.       double e1abs = fabs (e1);
  106.       double delta2 = e2 - e1;
  107.       double err2 = fabs (delta2);
  108.       double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
  109.       double delta3 = e1 - e0;
  110.       double err3 = fabs (delta3);
  111.       double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
  112.  
  113.       double e3, delta1, err1, tol1, ss;
  114.  
  115.       if (err2 <= tol2 && err3 <= tol3)
  116.     {
  117.       /* If e0, e1 and e2 are equal to within machine accuracy,
  118.          convergence is assumed.  */
  119.  
  120.       *result = res;
  121.       absolute = err2 + err3;
  122.       relative = 5 * GSL_DBL_EPSILON * fabs (res);
  123.       *abserr = GSL_MAX_DBL (absolute, relative);
  124.       return;
  125.     }
  126.  
  127.       e3 = epstab[n - 2 * i];
  128.       epstab[n - 2 * i] = e1;
  129.       delta1 = e1 - e3;
  130.       err1 = fabs (delta1);
  131.       tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
  132.  
  133.       /* If two elements are very close to each other, omit a part of
  134.          the table by adjusting the value of n */
  135.  
  136.       if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3)
  137.     {
  138.       n_final = 2 * i;
  139.       break;
  140.     }
  141.  
  142.       ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
  143.  
  144.       /* Test to detect irregular behaviour in the table, and
  145.          eventually omit a part of the table by adjusting the value of
  146.          n. */
  147.  
  148.       if (fabs (ss * e1) <= 0.0001)
  149.     {
  150.       n_final = 2 * i;
  151.       break;
  152.     }
  153.  
  154.       /* Compute a new element and eventually adjust the value of
  155.          result. */
  156.  
  157.       res = e1 + 1 / ss;
  158.       epstab[n - 2 * i] = res;
  159.  
  160.       {
  161.     const double error = err2 + fabs (res - e2) + err3;
  162.  
  163.     if (error <= *abserr)
  164.       {
  165.         *abserr = error;
  166.         *result = res;
  167.       }
  168.       }
  169.     }
  170.  
  171.   /* Shift the table */
  172.  
  173.   {
  174.     const size_t limexp = 50 - 1;
  175.  
  176.     if (n_final == limexp)
  177.       {
  178.     n_final = 2 * (limexp / 2);
  179.       }
  180.   }
  181.  
  182.   if (n_orig % 2 == 1)
  183.     {
  184.       for (i = 0; i <= newelm; i++)
  185.     {
  186.       epstab[1 + i * 2] = epstab[i * 2 + 3];
  187.     }
  188.     }
  189.   else
  190.     {
  191.       for (i = 0; i <= newelm; i++)
  192.     {
  193.       epstab[i * 2] = epstab[i * 2 + 2];
  194.     }
  195.     }
  196.  
  197.   if (n_orig != n_final)
  198.     {
  199.       for (i = 0; i <= n_final; i++)
  200.     {
  201.       epstab[i] = epstab[n_orig - n_final + i];
  202.     }
  203.     }
  204.  
  205.   table->n = n_final + 1;
  206.  
  207.   if (nres_orig < 3)
  208.     {
  209.       res3la[nres_orig] = *result;
  210.       *abserr = GSL_DBL_MAX;
  211.     }
  212.   else
  213.     {                /* Compute error estimate */
  214.       *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
  215.          + fabs (*result - res3la[0]));
  216.  
  217.       res3la[0] = res3la[1];
  218.       res3la[1] = res3la[2];
  219.       res3la[2] = *result;
  220.     }
  221.  
  222.   /* In QUADPACK the variable table->nres is incremented at the top of
  223.      qelg, so it increases on every call. This leads to the array
  224.      res3la being accessed when its elements are still undefined, so I
  225.      have moved the update to this point so that its value more
  226.      useful. */
  227.  
  228.   table->nres = nres_orig + 1;  
  229.  
  230.   *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));
  231.  
  232.   return;
  233. }
  234.